/* Copyright (C) Nial Peters 2014
 *
 * This file is part of pydoas.
 *
 * pydoas is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * pydoas is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with pydoas.  If not, see <http://www.gnu.org/licenses/>.
 */

#include"shift_and_squeeze.h"

#define FREE(x) if(x != NULL){free(x);}

///////////////////////////////////////////////////////////////////////////////

void
destroySNSDataObj(sns_data_t **data){

	int i; //counter
	sns_data_t *d = *data;

	if (d == NULL){
		return;
	}

	//free the arrays held in the struct
	for(i=0; i<d->num_refs; i++){

		if (d->accel_objs[i] != NULL){
			gsl_interp_accel_free(d->accel_objs[i]);
		}

		if (d->splines[i] != NULL){
			gsl_spline_free(d->splines[i]);
		}
	}

	FREE(d->accel_objs);
	FREE(d->splines);

	FREE(d->fitting_wavelengths);

	//now free the struct itself
	FREE(d);

	//finally set the pointer to NULL so that it cannot be used by mistake
	*data = NULL;
}

///////////////////////////////////////////////////////////////////////////////

sns_data_t *
createSNSDataObj(double **refs, int num_refs, double **ref_wavelengths,
				 int *ref_lengths, double *fitting_wavelengths,
				 int fit_window_length){

	int i; //counter
	sns_data_t *obj = NULL;

	if ((obj = malloc(sizeof(sns_data_t))) == NULL){
		ERROR(ERR_MEM, "Failed to allocate memory for shift and "
				       "squeeze data object.");
		return NULL;
	}

	//populate the struct with the values passed in
	obj->num_refs = num_refs;
	obj->fit_window_length = fit_window_length;

	//copy the fitting wavelengths array into the struct
	if ((obj->fitting_wavelengths = malloc(fit_window_length * sizeof(double))) == NULL){
		ERROR(ERR_MEM, "Failed to allocate memory to copy fitting wavelengths array.");
		destroySNSDataObj(&obj);
		return NULL;
	}

	memcpy(obj->fitting_wavelengths, fitting_wavelengths,
			fit_window_length * sizeof(double));

	//create the array for the accelerator objects
	if ((obj->accel_objs = malloc(num_refs * sizeof(gsl_interp_accel*)))==NULL){
		ERROR(ERR_MEM, "Failed to allocate memory for spline accelerator object array.");
		destroySNSDataObj(&obj);
		return NULL;
	}

	//create the array for the spline objects
	if ((obj->splines = malloc(num_refs * sizeof(gsl_spline*)))==NULL){
		ERROR(ERR_MEM, "Failed to allocate memory for spline object array.");
		destroySNSDataObj(&obj);
		return NULL;
	}

	//populate the splines and accelerator object arrays
	#pragma omp parallel for
	for (i=0; i<num_refs; i++){
		obj->accel_objs[i] = gsl_interp_accel_alloc();

		obj->splines[i] = gsl_spline_alloc(gsl_interp_cspline, ref_lengths[i]);

		if (obj->splines[i] != NULL){

			if (gsl_spline_init(obj->splines[i], ref_wavelengths[i], refs[i],
				ref_lengths[i]) != 0){

				//then we have failed to initialise the spline - set it to NULL
				gsl_spline_free(obj->splines[i]);
				obj->splines[i] = NULL;
			}
		}
	}

	//check that all the spline and accelerator objects were created successfully
	for (i=0; i<num_refs; i++){
		if ((obj->splines[i] == NULL) || (obj->accel_objs[i] == NULL)){
			ERROR(ERR_MISC, "Failed to create or initialise one or more spline objects.");
			destroySNSDataObj(&obj);
			return NULL;
		}
	}

	return obj;
}

///////////////////////////////////////////////////////////////////////////////

int
shiftAndSqueeze(const sns_data_t * sns_data, const double shift,
				const double squeeze, const int ref_num, double *output){

	int i;
	double lambda;

	//calculate the wavelength at the centre of the fitting window
	const double mid_point = 0.5 * (sns_data->fitting_wavelengths[0] + sns_data->fitting_wavelengths[sns_data->fit_window_length - 1]);
	const double offset = mid_point + shift;

	//check that the requested reference is valid
	if ((ref_num < 0) || (ref_num >= sns_data->num_refs)){
		ERROR(ERR_VAL, "Request for invalid reference '%d'. The shift and squeeze data object contains %d references.", ref_num, sns_data->num_refs);
		return 1;
	}

	//check that the requested shift and squeeze are within the interpolation limits
	lambda = (squeeze * (sns_data->fitting_wavelengths[0] - offset)) + mid_point;
	if (lambda < sns_data->splines[ref_num]->interp->xmin){
		ERROR(ERR_VAL, "The requested shift (%lf) and squeeze (%lf) values are outside of the interpolation limits for this reference and these fitting wavelengths.", shift, squeeze);
		return 1;
	}

	lambda = (squeeze * (sns_data->fitting_wavelengths[sns_data->fit_window_length - 1] - offset)) + mid_point;
	if (lambda > sns_data->splines[ref_num]->interp->xmax){
		ERROR(ERR_VAL, "The requested shift (%lf) and squeeze (%lf) values are outside of the interpolation limits for this reference and these fitting wavelengths.", shift, squeeze);
		return 1;
	}

	//compute the shifted and squeezed values within the fitting window
	for (i=0; i<sns_data->fit_window_length; i++){
		lambda = (squeeze * (sns_data->fitting_wavelengths[i] - offset)) + mid_point;
		output[i] = gsl_spline_eval(sns_data->splines[ref_num], lambda, sns_data->accel_objs[ref_num]);
	}

	return 0;

}

///////////////////////////////////////////////////////////////////////////////
